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Monte Carlo study of carrier-light coupling in terahertz quantum cascade 
lasers 
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We present a method for self-consistently including the optical cavity field into Monte Carlo-based carrier 
transport simulations. This approach allows for an analysis of the actual lasing operation in quantum cascade 
lasers, considering effects such as gain saturation and longitudinal mode competition. Simulation results for 
a terahertz quantum cascade laser are found to be consistent with experiment. 
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The development of innovative types of quantum cas- 
cade lasers (QCLs) and subsequent design optimization 
has gone hand in hand with detailed modelling, involv- 
ing more and more sophisticated simulation tools^— 
Especially in the terahertz regime, there is still plenty 
of room for improvement of the structures, for exam- 
ple in terms of output power, efficiency and temperature 
performance In this context, detailed carrier transport 
simulations have proven very useful, accounting for inter- 
and intrasubband processes alike and yielding not only 
level occupations, but also kinetic carrier distributions 
within each of the levels. Well-established approaches 
include the semiclassical ensemble Monte Carlo (EMC) 
method,— ~— and quantum transport simulations based on 
the density matris^ or non-equilibrium Green's functions 
formalismj^iili Notably, these simulations have up to now 
almost exclusively focused on the carrier transport, com- 
pletely neglecting the optical cavity field. An exception 
is a Monte-Carlo-based study of the coupled cavity dy- 
namics and electron transport, however only considering 
carrier interaction with photons and phonons. 1 ^ Carrier 
transport simulations allow for an analysis of the unsat- 
urated optical gain, indicating if the investigated QCL 
structure will lase at all under the assumed conditions. 
However, no statements about the actual lasing oper- 
ation, including the emitted optical power, the electric 
current and the carrier distributions, can be made. Evi- 
dently, an inclusion of the lasing action in the simulation 
would be desirable for many practical reasons, above all 
the device optimization with respect to the output power 
and wall-plug efficiency. Furthermore, such an analy- 
sis could provide insight into the carrier dynamics on 
a microscopic level, yielding an improved understand- 
ing of the light-matter interaction in such devices. In 
the following, we introduce an approach which allows 
for a straightforward implementation of the laser field 
into EMC simulations, without a significant increase of 
the numerical load. The presented scheme considers the 
coupled cavity field and carrier transport dynamics in a 
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completely self-consistent manner, accounting for effects 
such as gain saturation and longitudinal mode competi- 
tion. We present simulation results for a terahertz QCL, 
which are found to be consistent with experimental data. 

Induced optical transitions between an initial and a 
final level, in the following denoted by i and j, are com- 
monly described in terms of a spectral power gain coef- 
ficient gij (u>) and the population change associated with 
the induced emission and absorption events j» 
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Here, (w) > corresponds to gain, and g^ (w) < 
indicates loss. The constants Zq and h denote the 
impedance of free space and reduced Planck constant, 
respectively. V and uq are the volume and refractive in- 
dex of the gain medium. The line shape as a function of 
the angular optical frequency u> is given by 



Cij (w) = - 
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where 7^ and = (Ei — Ej) /K denote the optical 
linewidth and resonance frequency of the transition. E%j 
and Pi.j are the level eigenenergies and occupations, and 
dij is the corresponding transition matrix element. Fur- 
thermore, / denotes the optical intensity. 

In a semiclassical EMC simulation, the carrier trans- 
port is modeled by a stochastic evaluation of the inter- 
and intrasubband scattering events for a large ensemble 
of discrete particles. Each carrier is at a given time de- 
scribed by its quantum state |i„,k„) (which then has an 
occupation p = 1), where i n denotes the subband and 
k„ is the in-plane wave vector of the nth carrier, with 
n = 1, . . . ,N. In such simulations, the physical quanti- 
ties are typically computed from the corresponding en- 
semble averages. The total gain is obtained from Eq. 
(jlap by summing over all carriers n and all available final 
states \j, k„) with conserved in-plane wave vector k„, 
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For periodic structures like QCLs, it is appropriate to 
evaluate transitions from states within a single central 
period to the available final states (also including those 
in adjacent periods). The simulated volume is then given 
by V — LpN/n s , where L p and n s are the period length 
and the sheet doping density per period. 

In EMC simulations, all scattering processes are 
evaluated based on quantum mechanically calculated 
Boltzmann-like scattering rates. The contribution from 
optically induced transitions can be computed from Eq. 
(fib"]). To account for effects such as mode competition 
and multimode lasing, we have to sum over all relevant 
longitudinal modes, characterized by their frequencies 
oj m and intensities I m - For a carrier sitting in state \i, k), 
the transition rate to a state |j, k) with occupation prob- 
ability /jk thus becomes 

n^j = ^-^ \dij | 2 (1 - f jk ) T m£ij W • (4) 

In principle, in EMC simulations the photon dynamics 
can be evaluated in analogy to the carrier transport, by 
stochastic sampling of the discrete photon population in 
each relevant modeJ^. However, the statistical fluctua- 
tions associated with the rare photon emission events 
make it difficult to obtain sufficient accuracy^ Such 
problems are here avoided by using classical intensities 
I m rather than discrete photon populations. The inten- 
sity evolution over a time interval 5 can be described by 

I m (t + 6) = I m (t) exp ([T (oj m ) g (oj m ) - a (w m )] cS/n ) , 

(5) 

where cS/uq is the propagation distance of the light in the 
gain medium during that time. In Eq. ([5]), S has to be 
chosen short enough that g (uj m ) can be considered con- 
stant, r denotes the confinement factor, c is the vacuum 
speed of light, and a — a m + a w is an effective loss coef- 
ficient containing both the mirror and waveguide loss.— 
The description of outcoupling losses by a distributed co- 
efficient a m is very common^ and works especially well 
for moderate output coupling at the facets. Interference 
effects between adjacent modes and counterpropagating 
waves are not considered in our approach. Such effects 
can become relevant for the multimode dynamics and co- 
herent instabilities in modelocked lasers*^ 

Our self-consistent 3D EMC simulation includes all 
the essential scattering mechanisms, also accounting 
for Pauli's exclusion principled Carrier-carrier interac- 
tions are implemented based on the Born approxima- 
tion, taking into account screening in the full random 
phase approximation. 4 Also considered is scattering with 
impurities as well as acoustic and longitudinal-optical 
phonons, including nonequilibrium phonon effects. Inter- 
face roughness is implemented assuming typical values 
of 0.12 nm for the mean height and 10 nm for the cor- 
relation length.— We simulate four periods of the struc- 
ture, using periodic boundary conditions for the first and 
last period^ The optical cavity dynamics is considered 
by including the photon-induced scattering contribution, 




FIG. 1. (Color online) Temporal evolution of the outcoupled 
power P™' at different longitudinal mode frequencies f m . 
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FIG. 2. (Color online) Unsaturated and saturated power ma- 
terial gain coefficient vs frequency. For comparison, the cavity 
loss is also indicated. 

Eq. in the carrier transport simulation and updating 
g (uj m ) and I m for each considered mode after sufficiently 
short time intervals S, using Eqs. Q and ([5]). The ho- 
mogeneous linewidth in Eq. ^ is self-consistently cal- 
culated based on lifetime broadening. 7 Stimulated emis- 
sion and absorption do not affect the optical coherence; 
the corresponding scattering events are thus ignored in 
the linewidth calculation. Inhomogeneous gain broad- 
ening arises if transitions at different frequencies con- 
tribute to the gain spectrum, which can result in mul- 
timode operation. The subband energies and wave func- 
tions are obtained by solving the Schrodinger-Poisson 
(SP) equation^ The EMC and SP simulations are then 
performed iteratively until convergence is reached. 

To validate our simulation approach, we compare the 
results to experimental data for a 3.0 THz resonant 
phonon depopulation design, consisting of 178 periods. 14 
For this structure, comprehensive temperature resolved 
experimental data are available, along with detailed spec- 
ifications of the device, including the optical cavity. 
The cavity length and gain medium cross section are 
L = 1.22 mm and A = 178 x 0.0539 /xm x 23 fim, re- 
spectively, and the facet reflectivity is R — 0.85; further- 
more, no = 3.8, a w = 18.7 cm , a m — 1.3cm , and 
r = 0.93. 14 First, the structure is investigated for a lat- 
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FIG. 3. (Color online) Simulated and measured maximum 
outcoupled optical power vs temperature. The inset shows 
the simulated optical power vs the applied bias for different 
lattice temperatures Tl. 

tice temperature of TL = 50 K at a bias of 10.7kV/cm, 
where our simulation yields the maximum output power. 
We consider modes between 2.5 THz and 3.5 THz, with a 
frequency spacing of df = cj (2Luq) = 0.0324 THz. They 
are seeded with a small initial intensity (l m (t = 0) = 
313 W/cm 2 ). For moderate outcoupling, the single facet 
power of a mode at frequency f m — uj m / (27r) is obtained 
as P™* = ^I m A(l — R)/T, where the factor i arises 
because I m encompasses both the forward and backward 
propagating wave in the cavity. In Fig. [TJ the temporal 
evolution of P™* ^ s shown for 15 of the modes, centered 
around the gain maximum. After a time of 10 ns, the 
lasing operation has nearly reached steady state, yield- 
ing a stationary output power of 10.0 mW which is basi- 
cally accumulated in a single mode. This is in agreement 
with the experimentally observed predominant single- 
mode behavior of this structure.— The corresponding 
experimental value, measured at a heat sink tempera- 
ture T s = HKii and uncorrected for the collection ef- 
ficiency of the Winston cone used in the measurement 
setup, is 2.6mW.- With an estimated collection effi- 
ciency of around 30%ji£ this value corresponds to the 
simulation result. The unsaturated (solid line) and sat- 
urated (dashed line) spectral gain profiles are shown in 
Fig. [2J and the cavity loss is also indicated for compari- 
son (dotted line). As expected for stationary lasing, the 
peak saturated gain is clamped at the cavity loss. 

Multimode simulations, as presented in Fig. [TJ are 
computationally quite tedious, since long simulation 
times are required to reach steady state. Thus, in 
the following, we a priori assume single-mode opera- 
tion as experimentally observed,— considering in our 
simulation only the mode with the maximum gain. In 
Fig. [31 the maximum outcoupled optical power (nor- 
malized to its value at 50 K) is shown as a function of 
the lattice temperature Tl , comparing simulation results 
(solid line) to experimental temperature resolved mea- 
surements (dashed line)^^ 7 - The overall agreement is 



good, with a slight deviation between the simulated and 
measured maximum operating temperatures (175 K vs 
158 K), which we mainly attribute to uncertainties in the 
exact value of a w r^ The experimental structure lases best 
at around 13 V for all temperatures^ Also the simula- 
tion yields a temperature insensitive optimum bias which 
is located at around 10.7kV/cm (see the inset of Fig. 
[3]), corresponding to a voltage drop of 10.3 V across the 
active region of the investigated structure. This differ- 
ence between the theoretical and experimental value can 
largely be explained by additional parasitic voltage drops 
in the experimental structure, especially in the contacts.- 
We have also successfully tested our approach for various 
other designs, such as a 4.4 THz high power QCL. 19 

In conclusion, a method is presented to straightfor- 
wardly include the optical cavity field into self-consistent 
EMC carrier transport simulations, allowing for the anal- 
ysis of the actual lasing operation in QCLs. Effects like 
gain saturation, gain clamping and mode competition are 
naturally accounted for. Comparisons to experimental 
data confirm the validity of our approach. 
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